NDLab

Beta Version

The aim of this tutorial is to show how to query and extract data, plug them in Pyhton tools and custom code

To install NDLab:
go to the GitHub page and read the How to install section

In some examples this Notebook makes use of two external packages: pandas and plotly. You can of course use NDLab without them, they just make plotting easier.
To install plotly:

    $ pip install plotly
    $ pip install ipywidgets
    $ jupyter labextension install jupyterlab-plotly

To install pandas:

    $ pip install pandas

(You will have to restart jupyer lab)

This tutorial has three parts:

• 1^st - Work with data
  deals with extraction, analysis, plotting, and exporting of the data
• 2^nd - Write your own code
  deals with functions for writing code and building algorithms, like decay energy balance, decay chains, etc...
• 3^rd - Miscellanous examples
  to be used as templates to build your own library

Let's start importing NDLab:

# NuclearDataLab
import ndlab as nl

# this is needed only for using autocompletion when typing
from ndlaborm import * 

import pandas as pd
import io

1^st - Work with data

Just type the properties you want. As example, for the nuclide charge radius type NUCLIDE.CHARGE_RADIUS

Practice
Use the autofilling to see what is availabe: start typing 'NU' and hit the tab key to load NUCLIDE, then add '.', to see all the fields.
For the moment no code will be run, it is just to test the input method

Here the list of entites one can query
NUCLIDE, LEVEL, GAMMA, L_DECAY, CUM_FY, IND_FY
plus the set of decay radiations:
DR_ALPHA, DR_GAMMA, DR_BETAM, DR_BETAP, DR_ANTI_NU, DR_NU, DR_X, DR_ANNHIL, DR_AUGER, DR_CONV_EL, DR_DELAYED, DR_PHOTON_TOTAL

Practice
Try the autofilling with some of entities to see what they provide.
For the moment no code will be run, it is just to test the input method

Entities can have links to other entities
For example `GAMMA.START_LEVEL` gives access to the start level of the gamma, and if you type `GAMMA.START_LEVEL.ENERGY` you will refer to its energy

In the picture below the fields providing the links are in green - to nuclide, and red - to level.

Practice 1
Try to reach the angular momentum J of a level that decays emitting beta radiation (DR_BETA)

to see the solution, go to the cell at the very bottom

Practice
use print(description(NUCLIDE)), or any of the other names, to get the details the available fields and links

Note:
- (Q) tells that there are also the _UNC and _LIMIT fields. For example SP, SP_UNC, SP_LIMIT
- (q) tells that there is the _UNC field. For example ABUNDANCE, ABUNDANCE_UNC
- (S) tells that the fields is a string. For example NUCID, or JP

print(description(NUCLIDE))

NUCLIDE: properties of the nuclide in its ground state

ABUNDANCE 	*(q)* natural abundance, in mole fraction
ATOMIC_MASS 	*(Q)* atomic mass, in micro AMU
BETA_DECAY_EN 	*(Q)* energy available for beta decay [keV]
BINDING_EN 	*(Q)* binding energy per nucleon [keV]
CHARGE_RADIUS 	*(Q)* Root-mean-square of the nuclear charge radius, expressed in fm. 
ELEM_SYMBOL 	*(S)* symbol of the element
MASS_EXCESS 	*(Q)* mass excess [keV]
N 		number of neutrons
NUC_ID 		*(S)* identifier, e.g. 135XE
QA 		*(Q)* alpha decay Q energy [keV]
QBMN 		*(Q)* beta- + n decay energy [keV]
QEC 		*(Q)* electron capture Q-value [keV]
S2N 		*(Q)* 2-neutron separation energy [keV]
S2P 		*(Q)* 2-proton separation energy [keV]
SN 		*(Q)* neutron separation energy [keV]
SP 		*(Q)* proton separation energy [keV]
Z 		number of protons

1.1   How to build a data request ...

You will be using the grammar described above to get the data.
You need to specify what fields to retrieve and what filter, if any, applies.
For example to get all the energies of the γ-rays that start from a level with J^π = 2 ⁺
`fields = " GAMMA.ENERGY "`
`filter = " GAMMA.START_LEVEL.JP = '2+' "` The filter can be empty, but that in some cases can lead to a lengthy retrieval

Attention:
Enclose the fields and filter between double commas " "

To take advantage of the autofilling, first write everything without commas, then place them when you are ready.
This happens in a notebook, whilst on command line and other developing environments the autofilling works also within commas

... and get the data

The simplest way to retrieve data is to call the nl.csv_data(fields, filter) function, or the nl.json_data(fields, filter) function

Practice 2
use the fields and filter variable above to call the csv_data function (place it in a print to have a nicer display), then try the json_data function

to see the solution, go to the bottom cell

Some rules:

1) Only one entity
In a list of fields, you must use the same entity:
fields = "NUCLIDE.Z, GAMMA.ENERGY" is not allowed, use instead
fields = "GAMMA.NUC.Z, GAMMA.ENERGY"

2) Two dots max
Do not use more than 2 dots:
DR_BETA.DAUGHTER_FED_LEVEL.NUC.ABUNDANCE is not allowed: 3 dots.
Often you can do the same with 2 dots, in this case
DR_BETA.DAUGHTER.ABUNDANCE

3) Multiple fields
must be seperated by comma

4) Multiple conditions
Use the logical operators and , or to join conditions in the filter. Matematical operators are = > < %
In general, one can use SQL operators allowed by SQLite

Tips:
You can build functions with the fields:
fields = "GAMMA.MIXING_RATIO / GAMMA.ENERGY"

You can give alias for convenience
fields = "GAMMA.MIXING_RATIO / GAMMA.ENERGY as m"

A complex example

Let's consider this picture from K. S. Krane Phys. Rev. C 10, 1197 (1974)

Nuclide : Z even, N even

**-------------------** Start level: Jp = 2+     (2^ndoccurrence)
   ↓
   ↓        Gamma:   E2 + M1 multipolarity
   ↓
**-------------------** End level: Jp = 2+     (1^stoccurrence)


Filter criteria
- gammas from even-even nuclides
- starting from the second Jp = 2+ level
- ending at first Jp = 2+ level
- having E2+M1 multipolarity
Fields
- Δ = | mixing ratio / energy | * 835 , and A

How to specify conditions
Condition on text fields must be placed within commas ':
GAMMA.START_LEVEL.JP = '2+'

Condition on numeric fields do not have commas:
GAMMA.ENERGY > 2000

In case of doubt, use print(description(GAMMA)) to see the field type

<div>

Practice 3
Retrieve the data to reproduce the picture
The "fields" variable, and the first condition of the "filter", are already given.

- Remember the and to join conditions
- Use the autofilling and then add " at the beginning and at the end

to see the solution go to the cell at the very bottom

fields = "GAMMA.NUC.Z as z, GAMMA.NUC.N as n , abs( GAMMA.MIXING_RATIO / GAMMA.ENERGY / 0.835 * 1000) as m"

filter  = "( GAMMA.NUC.Z % 2 = 0 ) and ( GAMMA.NUC.N % 2 = 0 ) " # gamma from even-even nuclides
filter += " and GAMMA.START_LEVEL.JP = '2+' and GAMMA.START_LEVEL.JP_ORDER = 2 "# starts from Jp = 2+ , 2nd occurrence
#filter +=    # ends at     Jp = 2+ , 1st occurrence
#filter +=    # E2 + M1 multipolarity


#print(nl.csv_data(fields, filter)[0:250])

1.2   Plotting

So far only the NDLab package has been used, this is enough to retrieve data and export them in CSV or Json
What follows shows how to embed NDLab in Pandas and Plotly

Run the cell below which loads the packages

import numpy as np  
import pandas as pd 
import matplotlib.pyplot as plt 
import matplotlib.image as mpimg

# Plotly
import plotly.graph_objects as go
import plotly.express as px
   

The Y axis of the Krane's picture above needs Log(Δ). Using pandas, numpy, and matplotlib data are prepared and plotted.

With fields and filter of the above practice, the following cell loads the data and shows part of the dataframe

fields = " GAMMA.NUC.Z as z, GAMMA.NUC.N as n, abs( GAMMA.MIXING_RATIO / GAMMA.ENERGY / 0.835 * 1000 ) as y "   # z, n, mixing / energy

filter  = "(GAMMA.NUC.Z % 2 = 0 ) AND (GAMMA.NUC.N % 2 = 0 ) "   # even-even nuclides
filter += " and GAMMA.START_LEVEL.JP = '2+' and GAMMA.START_LEVEL.JP_ORDER = 2 "# starts from Jp = 2+ , 2nd occurrence
filter += " and GAMMA.END_LEVEL.JP = '2+' and GAMMA.END_LEVEL.JP_ORDER = 1 "   # ends at     Jp = 2+ , 1st occurrence
filter += " and ( GAMMA.MULTIPOLARITY = 'E2+M1' or  GAMMA.MULTIPOLARITY = 'M1+E2' )" # E2 + M1 multipolarity

# load the data into a dataframe
df = nl.pandas_df(fields, filter, pd)
df

	z	n	y
0	42	58	9.975355
1	44	56	5.386445
2	46	56	3.429497
3	46	64	12.527247
4	48	64	6.893959
...	...	...	...
129	40	50	0.266848
130	42	52	2.411730
131	42	54	0.732317
132	38	58	3.461286
133	40	56	0.226757

134 rows × 3 columns

Now the data are plotted by calculating Z + N for the x, and log(Δ) the the y

plt.scatter( df["z"] + df["n"],  np.log(df["y"])*10, color='red' )

<matplotlib.collections.PathCollection at 0x7fe7d910f730>

It is worth to superimpose the new data with the old ones

plt.subplots(figsize=(28, 10))
plt.scatter(   df["z"] + df["n"],  np.log(df["y"])*10, color='red')
plt.imshow(mpimg.imread("../docs/_static/krane_2.png"), extent=[38, 205, -63, 79])

<matplotlib.image.AxesImage at 0x7fe8016a4700>

The paper claims that there are minima at closed shells
On a 2D plot it is difficult to see where the closed shells in N or Z are placed
The cell below uses plotly to plot a 3D picture where nuclides with closed shells are in red:

• A second dataset is loaded with only the closed shells by adding a condition to the filter
• The first and second dataset are plotted togeter

Plotly graphs are interactive, can be rotated and zoomed

# condition for closed-shell 
filter2 = filter + " and ( GAMMA.NUC.Z  in ( 2,8,20,28,50,82,126 ) or GAMMA.NUC.N in ( 2,8,20,28,50,82,126 ) )"

# load data
df2 = nl.pandas_df(fields, filter2, pd)

# plot 3D
fig = go.Figure(data=[
    # all nuclides
    go.Scatter3d(
        x=df["z"], y=df["n"], z=np.log10(df.y)*10, 
        mode='markers', marker=dict(size=8, color='blue')
    )
 ,   
    # the magic ones in red
    go.Scatter3d(
        x=df2["z"], y=df2["n"], z=np.log10(df2.y)*10, 
        mode='markers',marker=dict(size=8,color="red")
    )
  
])

# set layout properties
fig.update_layout(margin=dict(l=0, r=0, b=30, t=0) ,width=800, height=800,    
                  scene = dict(xaxis_title='Z',yaxis_title='N',zaxis_title='Delta') 
                 ).show()

Once you have a dataframe, you can export into almost anything df.to_csv() df.to_json() df.to_hdf() df.to_excel() df.to_sql() df.to_xml()

Note: do not focus on the format, the format does not matter. Focus instead on the data model, the data model matters

Performance test

Get the energy distribution of the entire set of gamma transitions:
The following steps are performed with a couple lines of code

• 260 000 gamma lines in the db are taken,
• a condition on each line is applied,
• the results are ordered,
• the results are binned,
• an interactive plot is shown

fields =    "GAMMA.ENERGY as energy"
condition = "GAMMA.NUC.Z < 80 and  GAMMA.NUC.Z > 0 and GAMMA.ENERGY < 2000"

dfg =  nl.pandas_df(fields,condition,pd)

px.histogram( dfg,  x="energy", nbins=200).show()

Look at the binning above, is there anything relevant for clinical imaging of PET nuclides ?

1.3   Data analysis

The picture above relies on your brain to find patterns and outliers, whilst the cell below shows how to automatise the process.
The assumption is that the curve should be smooth, meaning that the first dervative should not have sudden jumps.
Comments in cell explain the steps

fields =    "GAMMA.ENERGY as energy"
condition = "GAMMA.ENERGY BETWEEN 0 AND 2000"

# load a pandas dataframe using the convenience funtion defined above
a =  nl.pandas_df(fields,condition,pd)

# bin the data (consult the web for further examples and explanations)
dft = pd.cut(a[a.columns[0]], bins=300).apply(lambda x: x.mid).value_counts(sort=False).to_frame()
dft = dft.rename(columns={"energy": "counts"})

# get the derivative
dft["delta"] = abs(dft.diff(periods = 1).counts )

# visualise the outlier (consult the web to check Plotly tools for data processing)
fig = px.box(dft, y=dft.delta).show()

# print the energy where the outlier occurs
print("Outlier energy ", dft.index[dft.index.get_loc(dft.delta.idxmax()) - 1 ])

Outlier energy  510.00550000000004

A few examples of data retrieval with filter to conclude part one of the tutorial

1.3.1   Plot log_ft values grouped by transition type

In the plot the transition type is not decoded: 1NU stands for 1^st non-unique, etc...

fields = "DR_BETA.TRANS_TYPE as t , DR_BETA.LOGFT as l"
filter = " DR_BETA.LOGFT != 0" # this is just to avoid empty values

df =  nl.pandas_df(fields, filter,pd)

# group by transition type
ans = [pd.DataFrame(y) for x, y in df.groupby(df.columns[0], as_index=False)]


# plot log_ft for each transition type 
fig = go.Figure()
for a in ans:
    tst = pd.cut(a['l'], bins=120).apply(lambda x: x.mid).value_counts(sort=False)
    fig.add_trace(go.Scatter(x=tst.index.values, y = tst, mode="markers",name=a.iloc[0, 0]))

fig.update_layout(width=800, height=400, xaxis_range=[3,15] ).show()

1.3.2   Plot log ft for allowed transitions form a parent level with Jp = '0+'

fields =  " DR_BETA.LOGFT as log_ft"
filter =  "( DR_BETA.PARENT_LEVEL.JP = '0+') and DR_BETA.TRANS_TYPE ='A'"

df = nl.pandas_df(fields, filter, pd)

px.histogram(df, x=df.log_ft, nbins=80).update_layout(width=350, height=350).show()

1.3.3   Extract level information for all even-even nuclei with the first excited state is NOT 2+.

With pandas it is easy to dump the data in almost any format. In this way one is then free to use them as input for other tools
Check the web for how to use the various writing - reading methods listed below

filter = ("LEVEL.NUC.N % 2 = 0 AND LEVEL.NUC.Z % 2 = 0 AND LEVEL.SEQNO = 1 AND LEVEL.JP != '2+' AND LEVEL.JP_METHOD = JP_STRONG ")
 
df = nl.pandas_df_nl(nl.levels(filter), pd)

df

	z	n	nucid	l_seqno	energy	energy_unc	energy_limit	half_life	half_life_unc	half_life_limit	...	jp_str	quadrupole_em	quadrupole_em_unc	quadrupole_em_limit	dipole_mm	dipole_mm_unc	dipole_mm_limit	questionable	configuration	isospin
0	36	64	100KR	1	309.00	10.0000	=	0.00	0.00	=	...	(2+)	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN
1	38	62	100SR	1	129.18	0.0009	=	3.91	0.16	=	...	(2+)	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN
2	38	64	102SR	1	126.00	0.2000	=	3.00	1.20	=	...	(2+)	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN
3	40	64	104ZR	1	139.30	0.0300	=	2.00	0.30	=	...	(2+)	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN
4	52	54	106TE	1	664.80	0.0300	=	0.00	0.00	=	...	(2+)	0.0	0.0	=	0.0	0.0	=	Y	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
155	74	118	192W	1	219.00	0.0000	=	0.00	0.00	=	...	[2+]	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN
156	60	100	160ND	1	65.20	0.0500	=	0.00	0.00	=	...	(2+)	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN
157	62	102	164SM	1	69.00	1.0000	=	0.00	0.00	=	...	(2+)	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN
158	92	124	216U	1	2240.00	42.0000	=	0.70	0.80	=	...	(8+)	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN
159	18	12	30AR	1	700.00	0.0000	=	0.00	0.00	=	...	(2+)	0.0	0.0	=	0.0	0.0	=	NaN	NaN	NaN

160 rows × 28 columns

2^nd - Functions and programming

Often one wants to write had hoc processing code that differs from filtering, grouping, binning, or any other manipulation. For example, calculate the energy balance of a decay.

Ndlab provides as a set of python classes and functions pushing the data extraction to the background. The functions accept the filter parameter (default no filter) to restrict the selection.

Note: Only the NDLab packages are needed in the following. One needs Pandas or Plotly just when using a dataframe or a plot

Warning: the filter can refer only to the entity retrieved: dr_gammas(" DR_GAMMA.ENERGY > 2000 ")

Tip: use the remove_doublers function in case the selection returns double entries :
delayed_n = nl.dr_delayeds( "DR_DELAYED.TYPE = DELAY_N ") parent_nucs = [ ld.parent for ld in delayed_n] parent_nucs = nl.remove_doublers(parent_nucs)
Whitout this call, a parent emitting many neutron will appear multiple times

The following is a list of hands-on examples. To see the full list of classes and function, please consult the guide

2.1   Decay energy balance

With NDLab's set of classes and functions, you do not need to take care of the data retrieval. See how you can calculate the energy balance of a decay.

filter = ""

# load a nuclide using its identifier
nuc = nl.nuclide("82SR")


# the decays() function loads an array with all the decays. Take the first one
decay = nuc.decays[0]

# all possible radiation types for the chosen decay
rads = [
        decay.xs(),          # X-ray
        decay.gammas(),      # Gamma
        decay.convels(),     # Conversion electron
        decay.augers(),      # Auger
        decay.alphas(),      # Alpha
        decay.betas_m(),     # B-
        decay.betas_p(),     # B+
        decay.anti_nus(),    # anti-neutrino
        decay.nus(),         # neutrino
        decay.annihil()      # annihilation
       ]

# measured energy for each radiation type
# note that for neutrini there are two contributions:
#      - one for the beta+ process, included in the first term below
#      - one for the electron capture process, included in the second
rads_en = [(sum([r[i].energy * r[i].intensity for i in range(len(r))])) for r in rads] + [(nu.energy_ec * nu.intensity_ec) for nu in decay.nus()]

# total measured energy
tot_en = sum(rads_en) / 100 

# intensities are per 100 decay of the parent, need to renormalise Q with branching ratio
q_br = decay.q_togs * decay.perc/100

# difference with deposited energy
delta = ((q_br - tot_en)/q_br)*100

print("Branching Ratio:          ", decay.perc/100 )
print("Energy accounted for [keV]" , tot_en )
print("Qb- * B.R. [keV]          " , q_br)
print("Delta [%]                  " , int(delta.value * 100) /100)
print("\nNotice the built-in uncertainty propagation")

Branching Ratio:            1.0+/-0
Energy accounted for [keV]  13.31+/-0.05
Qb- * B.R. [keV]            178+/-7
Delta [%]                   92.52

Notice the built-in uncertainty propagation

2.2  Ancestors - offsprings chain

Summary of decay chain, see the next cell for a more sophisticated example

# the daughters and parents of a nuclide are lists with the direct daughters and parents, respectively
# the list contains the Ground States of the daughters/parents

nuc = nl.nuclide("238U") 

print('Daughters and H-l [s]')
[print(d.nucid, d.gs.half_life_sec) for d in nuc.daughters]
    

print('\nParents and H-l [s]')
[print(p.nucid, p.gs.half_life_sec) for p in nuc.parents]


# the daughters_chain and parents_chain of a nuclide are lists with all daughters and parents, respectively
# the list contains the Ground States of the daughters/parents

print('\nDaughters chain \nNote that U-238 is listed as a daughter since it has an isomeric state that decays to the GS ')
dau = [n.gs for n in nl.nuclide("238U").daughters_chain]
for d in dau:
    print(d.nucid)


print('\nParents chain') 
par = [n.gs for n in nl.nuclide("205TL").parents_chain]
for p in par:
    print(p.nucid)
    

Daughters and H-l [s]
234TH  (2.0822+/-0.0026)e+06
238U  (1.4100+/-0.0019)e+17

Parents and H-l [s]
238PA  137+/-6
238U  (1.4100+/-0.0019)e+17
242PU  (1.183+/-0.006)e+13

Daughters chain 
Note that U-238 is listed as a daughter since it has an isomeric state that decays to the GS 
234TH
238U
234PA
234U
230TH
226RA
222RN
218PO
214PB
218AT
214BI
210TL
214PO
210PB
209PB
206HG
210BI
206TL
206PB
210PO
209BI
205TL
218RN

Parents chain
205HG
205PB
209BI
205AU
206PT
205PT
205IR
205BI
209PO
205PO
209AT
205AT
209RN
205RN
209FR
209RA
205FR
213TH
217U
209AC
213PA
213AC
217PA
213RA
217TH
221U
213FR
217AC
221PA
225NP
229AM
233BK
213RN
217RA
221TH
225U
229PU
233CM
237CF
241FM
209PB
213AT
223RA
210TL
209TL
213PO
223AC
223FR
227TH
227PA
231NP
231PU
235AM
239BK
239CF
243ES
243FM
247MD
223RN
227AC
223AT
223PO
223BI
224BI
227RA
231PA
231TH
227FR
231AC
235U
231RA
231FR
235PA
235NP
239PU
235TH
235PU
239AM
239CM
243CF
247FM
251NO
255RF
259SG
263HS
267DS
243BK
247ES
251MD
255LR
259DB
239NP
243CM
239U
243AM
239PA
243PU
247BK
244AM
243NP
247CM
243U
247AM
251CF
247PU
251BK
251ES
255FM
251CM
255ES
259MD
255CF
259NO
263RF
251FM
255MD
255NO
259RF
263SG
267HS
271DS
259LR
263DB
267BH
247CF
227RN
227AT
227PO
231U
211HG
214BI
214PB
218AT
215TL
214TL
218PO
216HG
215HG
214HG
218BI
222RN
218PB
222AT
226RA
226FR
226AC
230TH
226RN
226AT
226PO
230PA
230AC
234U
230RA
230FR
234PA
234NP
238PU
234TH
234AC
238U
234RA
238PA
242PU
238TH
242NP
242AM
246CM
242U
246AM
246BK
250CF
246PU
250CM
254CF
254ES
258MD
262LR
266DB
270BH
274MT
278RG
282NH
246CF
250ES
246ES
250FM
246FM
250MD
246MD
250NO
254RF
258SG
254LR
258DB
262BH
266MT
254NO
258RF
262SG
266HS
270DS
250BK
254FM
254MD
258LR
262DB
266BH
270MT
274RG
278NH
234PU
238AM
234AM
238CM
234CM
238CF
234BK
242CF
238BK
242ES
238NP
242CM
242BK
213BI
209HG
217AT
213PB
217PO
221FR
217BI
221RN
217PB
217TL
221AT
225RA
221PO
221BI
222BI
225FR
229TH
225RN
225AT
229RA
225PO
229FR
229RN
229AT
229AC
229PA
233U
229U
233NP
229NP
233PU
233AM
237CM
241CF
241ES
245FM
245MD
237AM
241BK
245ES
249MD
253LR
257DB
261BH
233PA
237PU
237NP
233TH
237U
241AM
237PA
241PU
237TH
241NP
245CM
241U
245AM
245BK
249CF
245PU
249BK
249CM
253ES
253CF
257FM
257ES
257MD
257NO
257LR
261RF
261DB
257RF
261SG
265HS
269DS
265SG
269HS
273DS
277CN
253FM
253MD
253NO
245CF
249ES
249FM
241CM
233AC
233RA
233FR
225AC
225TH
213TL
213HG
210AU
209AU
217RN
221RA
217FR
221AC
225PA

2.2.1   Detailed analysis

The function decay_chain(nucid, level, initial population), at the bottom of the next cell, produces a description of a decay chain.
The relevant entities are saved in a data structure for further analysis, see the final printout.

Note that in the Th-234 step, one needs to consider the Pa-234 metastable feeding
Indeed the Th-234 -> Pa-234 decay has two lines (#2 and #15), as well as the Pa-234 -> U-234 (#3 and #16)

This example shows again how NDLab allows to focus on the Physics of the issue, and not on libraries' formats, structure, conventions ...

# dictionary filled by the decay_chain() function to keep the relevant entities for further analysis
SAVE = []

# Keeps track of the decays already processed
DONE = []

# the H-l threshold for considering a metastable, in seconds
H_L_META = "60"

# do not consider dacays below thisthershold, in %
DECAY_PERC_TRESHOLD = 5

# filter any alpha, beta the decay radiation feeding a metastable : the seqno of the fed level must not be 0, and its H-l must be above the threshold
FILTER_META = "( DECAY_RAD.DAUGHTER_FED_LEVEL.HALF_LIFE_SEC > " + H_L_META +  " and DECAY_RAD.DAUGHTER_FED_LEVEL_SEQNO != 0 )"
# same as above but for gamma from IT decays
FILTER_META_GAMMA = "( DR_GAMMA.END_LEVEL.HALF_LIFE_SEC > " + H_L_META +  " and DR_GAMMA.END_LEVEL_SEQNO != 0 )"

# keeps track of nuclides already populated, to renormalise the production rate 
ALREADY_EXISTING = {}

# whether to print the decay balance
PRINT_BALANCE = True

# whether any radiation populates a metastable, consider only alpha, beta +/-, and isomeric transition
def get_rads_to_meta(d):
    
    if(d.code in [DECAY_A]):
            rads = d.alphas(FILTER_META)
    elif (d.code in [DECAY_Bm]):
            rads = d.betas_m(FILTER_META)
    elif (d.code in [DECAY_Bp]):
            rads = d.betas_p(FILTER_META)
    elif (d.code in [DECAY_IT]):
            rads = d.gammas(FILTER_META_GAMMA)
    else:
            print(' wrong decay ',d.code) # decay modes like delayed emission are not considered
            rads = []
            
    if(len(rads) > 1): print(rads[0].parent_nucid) # allow only 1 meta , warning if more 
    return rads

# the gamma feeding of a level from the above levels. For metastable, this intensity needs to be added to the direct feeding
def gamma_feed(decay, end_level):
    filter_gm = "DR_GAMMA.END_LEVEL_SEQNO = " + str( end_level)
    feeding = sum( g.intensity for g in decay.gammas(filter_gm))
    return feeding

# energy balance of the decay, to see whether the decay schema is complete
def energy_balance(decay):
     
    tot_en = decay.tot_measured_en()
    if(tot_en == 0): 
        return False
    print('\n', 'Decay',  decay.mode.name,'from', decay.nucid, 'level ',decay.l_seqno, 'to', decay.daughter_nucid , str(decay.perc)+' %')#, ' cumulative %',str((cum_intensity*100))+' %' )
    if( not PRINT_BALANCE): return True

    # intensities are per 100 decay of the parent, need to renormalise Q with branching ratio
    q_br = decay.q_togs * decay.perc/100

    # difference with deposited energy
    delta = ((q_br - tot_en)/q_br)*100
    
    print("Branching Ratio:          ", decay.perc/100 )
    print("Energy accounted for [keV]" , tot_en )
    print("Qb- * B.R. [keV]          " , q_br)
    print("Delta [%]                  " , int(delta.value * 100) /100)
    
    return True

# this performs the business
# it is called recursively from parent to daughter(s). A daughter can be processed more than once if there are many pathways to it
# nucid , l_seqno are from the  daughter of the previous decay, that now becomes the parent
# cum_intensity is the cumulative intensity, normalised to 100, with which the nuclide is produced:
def decay_chain(nucid, l_seqno, cum_intensity):
    
    # check if the cumulative intensity has a contribution from another decay already processed
    if nucid in ALREADY_EXISTING:
        cum_intensity = ALREADY_EXISTING[nucid]
   
    cum_intensity =  (cum_intensity/100)
    
    # get all decays
    nuc = nl.nuclide(nucid)
    decays = nuc.decays
   
    if( len(decays)==0) : return
    decays_todo = []
    
    # process only the decays from the given level
    for d in decays:
        if(d.perc > DECAY_PERC_TRESHOLD and d.l_seqno == l_seqno): decays_todo.append(d)
        
    for d in decays_todo:
        # skip if the decay was already processed
        
        if(d.pk in DONE or not energy_balance(d)): return    
        print('Cumulative population of the father:',str((cum_intensity*100))+' %' )
        # get the decay to a metastable
        rads_meta = get_rads_to_meta(d)
         
        # total intesitiy that populates the metastable : direct feeding plus gamma from above   
        tot_meta_feed = 0
        for r in rads_meta: 
            # no metastable populated, skip
            if(len(rads_meta) == 0): break
                
            tot_meta_feed = gamma_feed(d,r.fed_level.l_seqno) + r.intensity
            
            print('    metastable fed from parent level', r.parent_l_seqno, ' -->  to daughter level ',r.fed_level.l_seqno, ' gamma feeding: ' ,gamma_feed(d,r.fed_level.l_seqno) ,' direct feeding: ', r.intensity)
            
            # check if meta and the GS have the same daughter
            # in case, save the total BR for the daughter: when the daughter is processed it has the correct cumulative intensity
            dau_decays = nl.nuclide(d.daughter_nucid).decays   
            dau_dau = ''
            for dd in dau_decays:
                if(dd.l_seqno == 0): # decay from GS of the daughter
                    dau_dau = dd.daughter_nucid
                if(dd.l_seqno == r.fed_level.l_seqno and dd.daughter_nucid == dau_dau): # daughter of meta = daughter of GS
                    # simplify assuming that the meta has 2 decay modes, one IT to the GS, and one to the the same daughter of the GS
                    # example 234PA. More complex cases are not analysed here
                    ALREADY_EXISTING[dau_dau] = d.perc * cum_intensity
            
            # remember that the decay has been processed  
            DONE.append(d.pk)
            SAVE.append({'parent': nuc, 'level': l_seqno, 'fed_level':r.fed_level.l_seqno,'cumulative':(tot_meta_feed * cum_intensity), 'decay': d})
            # go to the next step in the chain using the decay from the meta
            decay_chain(d.daughter_nucid, r.fed_level.l_seqno, tot_meta_feed * cum_intensity)

        
        # for the GS take the decay BR minus what goes to a meta
        tot_gs_feed = d.perc - tot_meta_feed
        
        # remember that the decay has been processed  
        DONE.append(d.pk)
        SAVE.append({'parent': nuc, 'level': l_seqno, 'fed_level':0,'cumulative':(tot_gs_feed * cum_intensity), 'decay': d})
        
        # go to the next step in the chain using the decay from GS
        decay_chain(d.daughter_nucid, 0,tot_gs_feed * cum_intensity)    
        

###########################
# CALL THE FUNCTION 
# 238U is the initial nuclide, 
# 0 is its decaying level seqno (0 = Ground State)
# 100 is the initial normalisation of the number of nuclides

decay_chain('238U',0, 100)   
print('\nThe Dictionary has saved the relevant entities:\n')

cnt = 1
for s in SAVE:
    print('*'+str(cnt),s['parent'].nucid, s['decay'].daughter_nucid, s['level'],  s['fed_level'], s['cumulative'])
    cnt = cnt+1

 Decay A from 238U level  0 to 234TH  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV]  (4.20+/-0.18)e+03
Qb- * B.R. [keV]            4269.9+/-2.1
Delta [%]                   1.64
Cumulative population of the father: 100.0 %

 Decay B- from 234TH level  0 to 234PA  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV]  66.8+/-1.6
Qb- * B.R. [keV]            274.0+/-3.0
Delta [%]                   75.63
Cumulative population of the father:  100.0+/-0 %
    metastable fed from parent level 0  -->  to daughter level  2  gamma feeding:   4.23+/-0.28  direct feeding:   78.0+/-2.0

 Decay B- from 234PA level  2 to 234U  99.84+/-0.04 %
Branching Ratio:            0.9984+/-0.0004
Energy accounted for [keV] ? 1656.4+/-3.2
Qb- * B.R. [keV]            2264+/-4
Delta [%]                   26.84
Cumulative population of the father:  82.2+/-2.0 %

 Decay A from 234U level  0 to 230TH  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV] ? 4773+/-9
Qb- * B.R. [keV]            4857.5+/-0.7
Delta [%]                   1.73
Cumulative population of the father:  100.0+/-0 %

 Decay A from 230TH level  0 to 226RA  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV] ? 4679+/-19
Qb- * B.R. [keV]            4770.0+/-1.5
Delta [%]                   1.9
Cumulative population of the father:  100.0+/-0 %

 Decay A from 226RA level  0 to 222RN  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV]  4785+/-5
Qb- * B.R. [keV]            4870.70+/-0.25
Delta [%]                   1.76
Cumulative population of the father:  100.0+/-0 %

 Decay A from 222RN level  0 to 218PO  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV] ? 5485.1+/-0.5
Qb- * B.R. [keV]            5590.40+/-0.30
Delta [%]                   1.88
Cumulative population of the father:  100.0+/-0 %

 Decay A from 218PO level  0 to 214PB  99.980+/-0.020 %
Branching Ratio:            0.99980+/-0.00020
Energy accounted for [keV]  6001.34+/-0.16
Qb- * B.R. [keV]            6113.5+/-1.2
Delta [%]                   1.83
Cumulative population of the father:  100.0+/-0 %

 Decay B- from 214PB level  0 to 214BI  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV]  1016+/-10
Qb- * B.R. [keV]            1018+/-11
Delta [%]                   0.15
Cumulative population of the father:  99.980+/-0.020 %

 Decay B- from 214BI level  0 to 214PO  99.979+/-0.013 %
Branching Ratio:            0.99979+/-0.00013
Energy accounted for [keV] ? 3255+/-11
Qb- * B.R. [keV]            3268+/-11
Delta [%]                   0.41
Cumulative population of the father:  99.980+/-0.020 %

 Decay A from 214PO level  0 to 210PB  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV]  7686.82+/-0.09
Qb- * B.R. [keV]            7833.54+/-0.06
Delta [%]                   1.87
Cumulative population of the father:  99.959+/-0.024 %

 Decay B- from 210PB level  0 to 210BI  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV]  58.1+/-1.6
Qb- * B.R. [keV]            63.5+/-0.5
Delta [%]                   8.54
Cumulative population of the father:  99.959+/-0.024 %

 Decay B- from 210BI level  0 to 210PO  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV]  1161+/-4
Qb- * B.R. [keV]            1161.2+/-0.8
Delta [%]                   -0.01
Cumulative population of the father:  99.959+/-0.024 %

 Decay A from 210PO level  0 to 206PB  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV]  5304.39+/-0.07
Qb- * B.R. [keV]            5407.53+/-0.07
Delta [%]                   1.9
Cumulative population of the father:  99.959+/-0.024 %

 Decay B- from 234PA level  0 to 234U  100.0+/-0 %
Branching Ratio:            1.0+/-0
Energy accounted for [keV] ? (2.53+/-0.04)e+03
Qb- * B.R. [keV]            2194+/-4
Delta [%]                   -15.48
Cumulative population of the father:  17.8+/-2.0 %

The Dictionary has saved the relevant entities:

*1 238U 234TH 0 0  100.0+/-0
*2 234TH 234PA 0 2  82.2+/-2.0
*3 234PA 234U 2 0  82.1+/-2.0
*4 234U 230TH 0 0  100.0+/-0
*5 230TH 226RA 0 0  100.0+/-0
*6 226RA 222RN 0 0  100.0+/-0
*7 222RN 218PO 0 0  100.0+/-0
*8 218PO 214PB 0 0  99.980+/-0.020
*9 214PB 214BI 0 0  99.980+/-0.020
*10 214BI 214PO 0 0  99.959+/-0.024
*11 214PO 210PB 0 0  99.959+/-0.024
*12 210PB 210BI 0 0  99.959+/-0.024
*13 210BI 210PO 0 0  99.959+/-0.024
*14 210PO 206PB 0 0  99.959+/-0.024
*15 234TH 234PA 0 0  17.8+/-2.0
*16 234PA 234U 0 0  17.8+/-2.0

2.3   Level population fraction after a gamma cascade

Let's say that the level #15 of Xe-135 is populated by some reaction.

The following code gives the list of the levels populated by the gamma cascade, with energy and intensity

def cascade( my_level):
        global todo
        global levels
        global done
        
         
        for g in my_level.gammas() : 
            
            if(g.end_level.l_seqno != 0 and not (g.end_level.l_seqno in todo)):
                todo[g.end_level.l_seqno] = g.end_level
            
            if not g.end_level.l_seqno in result: 
                result[g.end_level.l_seqno] = result[my_level.l_seqno] * g.rel_photon_intens/100      
            else:
                result[g.end_level.l_seqno] += result[my_level.l_seqno] * g.rel_photon_intens/100         
        
        if not my_level.l_seqno in done: 
            done.append(my_level.l_seqno)
            
        todo = dict(sorted(todo.items(), reverse=True))  
        
        myrun = todo.copy()
       
        for r in myrun:
            if( not r in done):
                cascade(levels[r])

        return 



nuc_id = "135XE"
start_level = 15

levels = nl.nuclide(nuc_id).levels() 
result = {start_level : 1.0}
todo = {start_level : levels[start_level] }
done = []

cascade(levels[start_level])

# just printing
line = '|'.join(str(x).ljust(24) for x in [ "energy [keV] ", " population %", "level #"])
print(line +'\n')
for k in sorted(result):
    line = '|'.join(str(x).ljust(24) for x in [str(levels[k].energy ), str(result[k]*100.0),levels[k].pk])
    print(line)
  

energy [keV]            | population %           |level #                 

 0.0+/-0                | 100+/-8                |135XE-0                 
 288.455000+/-0.000015  | 0.0041+/-0.0007        |135XE-1                 
 526.551000+/-0.000013  | 34+/-5                 |135XE-2                 
 1131.512000+/-0.000011 | 7+/-5                  |135XE-3                 
 1260.416000+/-0.000013 | 0.137+/-0.024          |135XE-4                 
 1565.288000+/-0.000016 | 37+/-5                 |135XE-8                 
 1927.294000+/-0.000023 |100.0                   |135XE-15                

3^rd - Miscellaneous Examples

3.1  Photon total intensities

Nuclides often have many decay modes, either because of metastable states or because of a single energy state has more branchings

In gamma spectroscopy it is useful to have the total intensity of an energy line in the decay of a parent, regardless of the decay mode

The cell below selects Eu-152 and calls the dr_photon_total function which performs the task. Note that the function is called with a filter which cuts intensityes less than 2%

df = nl.pandas_df_nl(nl.nuclide("152EU").decays[0].dr_photon_tot("DR_PHOTON_TOTAL.INTENSITY > 2"), pd)

go.Figure().add_trace(go.Scatter(x=df['energy'], y = df['intensity'], mode = "markers")).update_layout(xaxis_title='Energy [keV]',
yaxis_title='Intensity [%]').show()

df

	parent_nucid	parent_l_seqno	energy	energy_unc	energy_limit	intensity	intensity_unc	intensity_limit	type	count
0	152EU	0	6.3540	0.0000	=	14.000	0.500	=	X	1.0
1	152EU	0	39.5220	0.0000	=	20.870	0.200	=	X	1.0
2	152EU	0	40.1170	0.0000	=	37.800	0.300	=	X	1.0
3	152EU	0	45.9980	0.0000	=	14.850	0.200	=	X	1.0
4	152EU	0	121.7817	0.0003	=	28.530	0.160	=	G	1.0
5	152EU	0	244.6974	0.0008	=	7.550	0.040	=	G	1.0
6	152EU	0	344.2785	0.0012	=	26.590	0.200	=	G	1.0
7	152EU	0	411.1165	0.0012	=	2.237	0.013	=	G	1.0
8	152EU	0	443.9606	0.0016	=	2.827	0.014	=	G	1.0
9	152EU	0	778.9045	0.0024	=	12.930	0.080	=	G	1.0
10	152EU	0	867.3800	0.0300	=	4.230	0.030	=	G	1.0
11	152EU	0	964.0570	0.0050	=	14.510	0.070	=	G	1.0
12	152EU	0	1085.8370	0.0100	=	10.110	0.050	=	G	1.0
13	152EU	0	1112.0760	0.0030	=	13.670	0.080	=	G	1.0
14	152EU	0	1408.0130	0.0030	=	20.870	0.090	=	G	1.0

3.2   Ground states Half-lives as function of N and Z

See that one can refer to a column in a dataframe by its position, in this way one has a template to plot 3d

# convenience function to plot 3D
def plot3d(fields, condition):
    df =  nl.pandas_df(fields,condition, pd)

    fig = go.Figure(data=[
    
    go.Scatter3d(
        x=df[df.columns[0]], y = df[df.columns[1]], z = df[df.columns[2]],
        mode='markers', marker=dict(size=8, color='blue') )

    ])

    fig.update_layout(height=1000 ,scene=dict(zaxis=dict( type='log'))).show()
    return df

fields =    "LEVEL.NUC.Z , LEVEL.NUC.N , LEVEL.HALF_LIFE_SEC"
condition = "LEVEL.ENERGY = 0"

plot3d(fields, condition)    

	z	n	half_life_sec
0	47	53	1.206000e+02
1	48	52	4.910000e+01
2	49	51	5.650000e+00
3	36	64	7.000000e-03
4	42	58	2.212141e+26
...	...	...	...
3349	92	122	5.200000e-04
3350	78	87	2.600000e-04
3351	80	90	8.000000e-05
3352	9	4	4.517205e-22
3353	17	35	NaN

3354 rows × 3 columns

3.3   Half-life of nuclides with 10 < Z < 100, emitting delayed neutrons after beta- decay

fields = "DR_DELAYED.PARENT.Z , DR_DELAYED.PARENT.N , DR_DELAYED.PARENT_LEVEL.HALF_LIFE_SEC"
condition = "DR_DELAYED.TYPE = 'DN' and DR_DELAYED.PARENT.Z > 10 and DR_DELAYED.PARENT.Z < 100"

plot3d(fields, condition)

	z	n	half_life_sec
0	37	63	0.0520
1	19	34	0.0300
2	49	84	0.1650
3	51	84	1.6790
4	33	52	2.0210
5	19	30	1.2600
6	29	50	0.2410
7	31	52	0.3081
8	19	31	0.4720
9	12	21	0.0905
10	19	33	0.1100
11	19	32	0.3650
12	37	56	5.8400
13	37	59	0.2030
14	37	58	0.3777
15	37	57	2.7020
16	47	76	0.2990
17	53	86	2.2800
18	53	84	24.5000
19	53	85	6.2600
20	35	53	16.3400
21	17	29	0.2320
22	35	55	1.9100
23	11	19	0.0480
24	35	54	4.3570
25	35	52	55.6800
26	11	18	0.0441
27	11	21	0.0132
28	35	57	0.3140
29	37	60	0.1691
30	37	62	0.0540

3.4   Half-life of nuclides in the decay chain of Am-241

For each nuclide, the properties daughters_chain and parents_chain contain the offsprings and the ancestors, respectively

The example shows the seamless interplay between sets of NDLab classes and pandas dataframes

# get the offsprings as NDLab class instances
dau = [n.gs for n in nl.nuclide("241AM").daughters_chain]

# dump in a dataframe
df = nl.pandas_df_nl(dau,pd)

# Plot. Hover on a point to see more info
px.scatter(df, 
           x= df.index, 
           y= df.half_life_sec,
           hover_data=["nucid"],
           labels={ "index": "Daughter # ", "y":"Half-life [s]"},
           log_y=True).show()

3.5  Level energy vs J

The filter is on Z = 20, and for J having a unique assignment

fields = "LEVEL.J as j, LEVEL.ENERGY as e"
filter =  "LEVEL.NUC.Z = 20 and LEVEL.JP_METHOD = RIPL_J_UNIQUE"

df = nl.pandas_df(fields, filter, pd).sort_values(by=['j'])

go.Figure().add_trace(go.Scatter(x=df['j'], y = df['e'], mode = "markers")).update_layout(xaxis_title='J',
yaxis_title='Energy').show()

3.6  Gamma energy vs J of the start level

fields = "GAMMA.START_LEVEL.J as j, GAMMA.ENERGY as e"
filter = " GAMMA.NUC_ID = '16O' " 

df = nl.pandas_df(fields, filter, pd).sort_values(by=['j'])

go.Figure().add_trace(go.Scatter(x=df['j'], y = df['e'], mode = "markers")).update_layout(xaxis_title='J',
yaxis_title='Energy').show()

3.7  Schmidt diagram for Z-odd, N-even nuclides

fields =  "LEVEL.J as j, LEVEL.DIPOLE_MM as m"
filter =  "LEVEL.NUC.Z % 2 = 1 and LEVEL.NUC.N % 2 = 0 and  LEVEL.JP_METHOD = RIPL_J_UNIQUE"

df = nl.pandas_df(fields, filter, pd).sort_values(by=['j'])

go.Figure().add_trace(go.Scatter(x=df['j'], y = df['m'], mode = "markers")).update_layout(xaxis_title='J',
yaxis_title='Magn. dipole [mu N]').show()

3.8  Linear relation between log(T1/2) and sqrt(Qval) for α decay of even-even nuclides

fields = "LOG(L_DECAY.LEVEL.HALF_LIFE_SEC) as h , 1/ SQRT(L_DECAY.NUC.QA + L_DECAY.LEVEL.ENERGY) as e "
filter =  "L_DECAY.NUC.Z = 92  and  L_DECAY.NUC.N % 2 = 0 and L_DECAY.MODE = DECAY_A "

df = nl.pandas_df(fields, filter, pd).sort_values(by=['e'])

go.Figure().add_trace(go.Scatter(x=df['e'], y = df['h'], mode = "markers")).update_layout(xaxis_title='E(keV)<sup>-1/2</sup>',
yaxis_title='log(T<sub>1/2</sub>)').show()

4^th - Solutions to the exercises

# practice 1
DR_BETA.PARENT_LEVEL.JP

# practice 2
fields = " GAMMA.ENERGY "
filter = " GAMMA.START_LEVEL.JP = '2+' "
print(nl.csv_data(fields, filter))
print(nl.json_data(fields, filter))

# practice 3
fields = "GAMMA.NUC.Z as z, GAMMA.NUC.N as n , abs( GAMMA.MIXING_RATIO / GAMMA.ENERGY / 0.835 * 1000) as b"

filter  = "( GAMMA.NUC.Z % 2 = 0 ) and ( GAMMA.NUC.N % 2 = 0) "                     # even-even nuclides
filter +=  " and GAMMA.START_LEVEL.JP = '2+' and GAMMA.START_LEVEL.JP_ORDER = 2 "   # starts from Jp=2+ ,2nd occurrence
filter +=  " and GAMMA.END_LEVEL.JP = '2+'   and GAMMA.END_LEVEL.JP_ORDER = 1 "     # ends at Jp = 2+ ,1st occurrence
filter +=  " and (GAMMA.MULTIPOLARITY = 'E2+M1'or GAMMA.MULTIPOLARITY = 'M1+E2') "  # E2 + M1 multipolarity

nl.csv_data(fields, filter)